#
# Plots.R
#
# Plots for
# Hill, Seth J., Daniel J. Hopkins, and Gregory A. Huber. 2021. 
# "Not by Turnout Alone: Measuring the Sources of Electoral Change, 2012 to 2016." Science Advances.
#

rm(list=ls())
if (.Platform$OS == "windows") {
  # Set working directory in location of script.
  .doit <- function() { # only works with R.exe; trap errors if using Rscript.exe
  f_f <- lapply(sys.frames(),function(x) x$ofile);f_f <- Filter(Negate(is.null),f_f) ; PTH <- dirname(f_f[[length(f_f)]]); setwd(PTH) ; rm(PTH,f_f)
  }
  try(.doit(),silent=T)
}

library(bit64)
library(data.table)
library(haven)
library(RColorBrewer)
palette(brewer.pal(n=9,name="Set1")[c(1:5,7:9)]) # reset default colors

DT = data.table(read_dta("PooledPrecinctAnalysisData.dta"))

#
# Recodes and functions.
#

# Compositional change.
# In counts.
DT[,GOPTURNOUTNUM12 := P12_V12_GOP + PBOTH_V12_GOP + PBOTH_VBOTH_GOP]
DT[,GOPTURNOUTNUM16 := P16_V16_GOP + PBOTH_V16_GOP + PBOTH_VBOTH_GOP]
DT[,DEMTURNOUTNUM12 := P12_V12_DEM + PBOTH_V12_DEM + PBOTH_VBOTH_DEM]
DT[,DEMTURNOUTNUM16 := P16_V16_DEM + PBOTH_V16_DEM + PBOTH_VBOTH_DEM]
DT[,comp.change := GOPTURNOUTNUM16 - GOPTURNOUTNUM12 - DEMTURNOUTNUM16 + DEMTURNOUTNUM12]
# In shares.
DT[,comp.change.shares := GOPTURNOUTPCT16 - GOPTURNOUTPCT12 - DEMTURNOUTPCT16 + DEMTURNOUTPCT12]
# Vote change in counts.
DT[,vote.change := GOP_Voters_2016_Prct - GOP_Voters_2012_Prct - Dem_Voters_2016_Prct + Dem_Voters_2012_Prct]
# Vote change in shares.
DT[,share.change := (GOPPCT16- DEMPCT16) - (GOPPCT12 - DEMPCT12)]

# Precinct stability := intersection/union.
DT[,union.electorate := P12_V12_DEM + P16_V16_DEM + PBOTH_V12_DEM + PBOTH_V16_DEM + PBOTH_VBOTH_DEM + P12_V12_GOP + P16_V16_GOP + PBOTH_V12_GOP + PBOTH_V16_GOP + PBOTH_VBOTH_GOP + P12_V12_OTH + P16_V16_OTH + PBOTH_V12_OTH + PBOTH_V16_OTH + PBOTH_VBOTH_OTH] 
DT[,intersection.electorate := PBOTH_VBOTH_DEM + PBOTH_VBOTH_GOP + PBOTH_VBOTH_OTH] 
DT[,stability := intersection.electorate/union.electorate]
print(DT[,summary(stability)])

.inRange <- function(x,probs=c(.025,.975),max.size=4500) {
  # Return boolean if x inside the quantiles in arg probs.
  flag = x >= quantile(x,probs[1]) & x <= quantile(x,probs[2])
  if (sum(flag) > max.size) {
    # Limit to max.size cases.
    samp = sample(which(flag),size=max.size,replace=F)
    flag = flag & seq(1,length(flag)) %in% samp
  }
  return(flag)
}

#
# %%%%%%%%%%%%%%%%%%%%%%%%%%
# Plot vote change on stability change.
# %%%%%%%%%%%%%%%%%%%%%%%%%%
#
.stabPlot <- function(DT,smth=T,span=1/4,y.var="vote.change") {
  # Make the plot.
  DT[,plot(x=stability,y=get(y.var),type='n',ann=F,axes=F,xlim=c(0,1))]
  axis(1);axis(2,las=2)
  grid(lwd=2)
  abline(h=0,col='black',lwd=3)
  DT[,points(x=stability,y=get(y.var),pch=19,cex=.7,col=rgb(0,0,0,alpha=.2))]
  title(xlab="Proportion of all voters voting twice",ylab="Electoral change toward GOP")
  if (smth) { # Add loess-smooth
    # Estimate smoother on further interior 95 pct.
    s = DT[.inRange(stability),loess.smooth(x=stability,y=get(y.var),span=span)]
    # Plot smoother.
    lines(x=s$x,y=s$y,lwd=5,col=1)
  }
}

# Counts.
pdf("ElectoralChangeByStability.pdf",height=6,width=6)
par(mar=c(4.1,4.1,1.6,0.6))
# Limit to interior 95 pct of data.
states = data.table(st=c("FL","GA","MI","NV","OH","PA"),st.nm=c("Florida","Georgia","Michigan","Nevada","Ohio","Pennsylvania"))
for (i in 1:nrow(states)) {
  states[ i, .stabPlot(DT[state == st,][.inRange(vote.change),]) ]
  states[ i, title(main=st.nm) ]
}
dev.off()

# Shares.
pdf("ElectoralChangeByStability-Shares.pdf",height=6,width=6)
par(mar=c(4.1,4.1,1.6,0.6))
# Limit to interior 95 pct of data.
states = data.table(st=c("FL","GA","MI","NV","OH","PA"),st.nm=c("Florida","Georgia","Michigan","Nevada","Ohio","Pennsylvania"))
for (i in 1:nrow(states)) {
  states[ i, .stabPlot(DT[state == st,][.inRange(share.change),],y.var="share.change") ]
  states[ i, title(main=st.nm) ]
}
dev.off()


#
# %%%%%%%%%%%%%%%%%%%%%%%%%%
# Plot vote change on compositional change.
# %%%%%%%%%%%%%%%%%%%%%%%%%%
#
.compPlot <- function(DT,smth=T,span=1/4,x.var="comp.change",y.var="vote.change") {
  # Make the plot.
  #lims = DT[,range(c(get(x.var),get(y.var)),na.rm=T)] # aspect=1
  lims = NULL
  DT[,plot(x=get(x.var),y=get(y.var),type='n',ann=F,axes=F,ylim=lims,xlim=lims)]
  axis(1);axis(2,las=2)
  grid(lwd=2)
  abline(v=0,col='black',lwd=3)
  abline(h=0,col='black',lwd=3)
  DT[,points(x=get(x.var),y=get(y.var),pch=19,cex=.7,col=rgb(0,0,0,alpha=.2))]
  title(xlab="Compositional change toward GOP",ylab="Electoral change toward GOP")
  if (smth) { # Add loess-smooth
    # Estimate smoother on further interior 95 pct.
    s = DT[.inRange(get(x.var)),loess.smooth(x=get(x.var),y=get(y.var),span=span)]
    # Plot smoother.
    lines(x=s$x,y=s$y,lwd=5,col=1)
  }
}

# Counts.
pdf("ElectoralChangeByPartisanCompositionChange.pdf",height=6,width=6)
par(mar=c(4.1,4.1,1.6,0.6))
# Limit to interior 95 pct of data.
states = data.table(st=c("FL","GA","MI","NV","OH","PA"),st.nm=c("Florida","Georgia","Michigan","Nevada","Ohio","Pennsylvania"))
for (i in 1:nrow(states)) {
  states[ i, .compPlot(DT[state == st,][.inRange(comp.change) & .inRange(vote.change),]) ]
  states[ i, title(main=st.nm) ]
}
dev.off()

# Shares.
pdf("ElectoralChangeByPartisanCompositionChange-Shares.pdf",height=6,width=6)
par(mar=c(4.1,4.1,1.6,0.6))
# Limit to interior 95 pct of data.
states = data.table(st=c("FL","GA","MI","NV","OH","PA"),st.nm=c("Florida","Georgia","Michigan","Nevada","Ohio","Pennsylvania"))
for (i in 1:nrow(states)) {
  states[ i, .compPlot(DT[state == st,][.inRange(comp.change.shares) & .inRange(share.change),],x.var="comp.change.shares",y.var="share.change") ]
  states[ i, title(main=st.nm) ]
}
dev.off()



